####################################### MULTIPLE EXPOSURE ##################################################


#### Opening the data ####

# You need to put your address
aids<-read.table("C:/Users/laupaz/Desktop/Teaching & visiting/Causal inference (Milano Bicocca) Dicembre 2018/Lab/multiple.txt",header=TRUE)
# Visualizing the first 10 rows
aids[1:10,]

#### Ignoring L1 ####
# Standard adjustment in R fitting generalized linear models with glm

summary(Standard<-glm(formula=Y~A0+A1+L0,family=binomial(link = "logit"),data=aids))

#### Considering L1 ####

# Sequential adjustment at t = 0 in R
summary(Sequentialt0<-glm(formula=Y~A0+L0,family=binomial,data=aids))

summary(Sequentialt1<-glm(formula=Y~A1+L0+A0+L1,family=binomial,data=aids))


#### Fitting of MSMs ####

# STEP 1 
# Fit a regression model for the exposure at each time point, given the observed past up to that time point
fit0=glm(formula=A0~L0,family=binomial,data=aids)
fit1=glm(formula=A1~L0+A0+L1,family=binomial,data=aids)

# STEP 2
# For each subject, use the fitted exposure model to estimate a subject-specific weight W = 1/ Pr(A0|L0)*Pr(A1|L0;A0;L1)

pred0=predict(object=fit0,type="respons")
pred1=predict(object=fit1,type="respons")
w0=1/(aids$A0*pred0+(1-aids$A0)*(1-pred0))
w1=1/(aids$A1*pred1+(1-aids$A1)*(1-pred1))
w=w0*w1

# STEP 3
# Fit the MSM using weighted regression, as if it would have been a model for Pr(Y = 1|A0;A1),
# summary(glm(formula=Y~A0+A1,family=binomial,data=aids,weights=w))
# glm warns if it detects that the no. of trials or successes is non-integral, 
# but it goes ahead and fits the model anyway. If you want to suppress the warning (and you're sure it's not a problem), 
# use family=quasibinomial instead

msm<-glm(formula=Y~A0+A1,family=quasibinomial,data=aids,weights=w)
summary(msm)
#  Use this more elaborate model to estimate the joint effect 
# of A0 and A1 (Beta0+Beta1) as los odds ratio
Joint_effect<-coef(msm)["A0"]+coef(msm)["A1"] # 2.345


#  (b) Add an interaction term between a0 and a1 in the MSM. 
msm1<-glm(formula=Y~A0+A1+A0*A1,family=quasibinomial,data=aids,weights=w)
summary(msm1)
# (b) Use this more elaborate model to estimate the joint effect 
# of A0 and A1 (Beta0+Beta1+Beta2) as log odds ratio
Joint_effect<-coef(msm1)["A0"]+coef(msm1)["A1"]+coef(msm1)["A0:A1"] # 1.5322


# Sandwich estimator to have robust standard errors considering the weigthing system
# Installing packages
library("sandwich", lib.loc="~/R/R-3.3.0/library")
library("lmtest", lib.loc="~/R/R-3.3.0/library")
model<-glm(Y~A0+A1,family=quasibinomial,data=aids,weights=w)
summary(model)
coeftest(model, vcov = sandwich) 

